---
title: "ACE"
author: "Edmund Kelly"
date: "2023-04-17"
output: html_document
---

```{r Packages, echo = FALSE} 
require(OpenMx)
library(stringr)
require(umx)
require(haven)
require(broom)
require(tidyverse)
require(stargazer)
require(eeptools)
require(ggplot2)
require(eeptools)
set.seed(1234)
```

```{r Pre-Processing}

#Loading data
STR <- data.frame(read_dta("1 - Data/cohort_combined.dta"))
Salty = data.frame(read_dta("1 - Data/slty_map.dta"))

combined = merge(Salty, STR, by.x = "LopNr", by.y = "LopNr")

#Removing missing values and incomplete pairs
combined[combined == 99] <- NA
combined[combined == 999] <- NA

#Trust

combined$oldtrust = ifelse(combined$ATTITYD31 == 9, NA, combined$ATTITYD31-1) 

combined$trust = (3-combined$oldtrust)/3

combined$oldtrust = 3-combined$oldtrust

combined$trust_scaled = as.numeric(scale(combined$trust))

#Behavioral

combined <- combined %>% mutate(
                          contactpol = 
                          ifelse(is.na(ATTITYD38_1) == F, yes = 1, no = 0), 
                          contactpub = 
                          ifelse(is.na(ATTITYD38_2) == F, yes = 1, no = 0), 
                          protest = 
                          ifelse(is.na(ATTITYD38_3) == F, yes = 1, no = 0), 
                          boycott = 
                          ifelse(is.na(ATTITYD38_4) == F, yes = 1, no = 0), 
                          finance = 
                          ifelse(is.na(ATTITYD38_5) == F, yes = 1, no = 0), 
                          petition = 
                          ifelse(is.na(ATTITYD38_6) == F, yes = 1, no = 0), 
                          behavioral = (contactpol + contactpub + protest + boycott + finance + petition), 
                          interpersonal_1 = ATTITYD2, 
                          interpersonal_2 = ATTITYD3, 
                          interpersonal_scaled = as.numeric(scale(interpersonal_1 + interpersonal_2)))

#Turnout
Turnout_1970 = data.frame(read_dta("1 - Data/Lev_valdelt_1970.dta"))
Turnout_1994f = data.frame(read_dta("1 - Data/Lev_valdelt_1994f.dta"))
Turnout_1994rkl = data.frame(read_dta("1 - Data/Lev_valdelt_1994rkl.dta"))
Turnout_2010 = data.frame(read_dta("1 - Data/vd10.dta"))
Turnout_2009 = data.frame(read_dta("1 - Data/vd09.dta"))
Turnout_2018 = data.frame(read_dta("1 - Data/Lev_vd18.dta"))
Turnout_2019_EU = data.frame(read_dta("1 - Data/EU2019.dta"))

Turnout_1994 = merge(Turnout_1994f, Turnout_1994rkl, by.x = "LopNr", by.y = "LopNr", all = TRUE)
Turnout_20th = merge(Turnout_1994, Turnout_1970, by.x = "LopNr", by.y = "LopNr", all = TRUE) %>% 
  rename("Ref_1994" = "f", 
         "Ref_1994_Abroad" = "utlandsrost.x",
         "Abroad_1994" = "utlandsrost.y",
         "Nat_1994" = "r.x", 
         "Reg_1994" = "l.x", 
         "Mun_1994" = "k.x", 
         "Nat_1970" = "r.y", 
         "Reg_1970" = "l.y", 
         "Mun_1970" = "k.y", 
         "Legal_1970" = "omyndig", 
         "Foreign_1970" = "utlmedb")

Turnout_0s = merge(Turnout_2009, Turnout_2010, by.x = "LopNr", by.y = "LopNr", all= TRUE) %>% 
  rename("EU_2009" = "e",  
         "Nat_2010" = "r", 
         "Reg_2010" = "l", 
         "Mun_2010" = "k", 
         "Eligible_2010" = "Rostratt.x")

Turnout_10s = merge(Turnout_2018, Turnout_2019_EU, by.x = "LopNr", by.y = "LopNr", all = TRUE) %>% 
  rename("EU_2019" = "Eurost", 
         "Nat_2018" = "rrost", 
         "Reg_2018" = "lrost", 
         "Mun_2018" = "krost", 
         "Eligible_2018" = "Rostratt.x")

Turnout_most = merge(Turnout_20th, Turnout_0s, by.x = "LopNr", by.y. = "LopNr", all = TRUE)

Turnout_all = merge(Turnout_most, Turnout_10s, by.x = "LopNr", by.y = "LopNr", all = TRUE) %>% 
  select(LopNr, EU_2019, Nat_2018, Reg_2018, Mun_2018, Eligible_2018, Eligible_2010, Mun_2010, Reg_2010, Nat_2010, EU_2009, Mun_1970, Reg_1970, Nat_1970, Mun_1994, Reg_1994, Nat_1994, Ref_1994_Abroad, Ref_1994, Abroad_1994, 
         Legal_1970, Foreign_1970)

combined = merge(Turnout_all, combined, by.x = "LopNr", by.y = "LopNr", all.y = TRUE) 

combined = combined %>% 
  mutate(EU_2009 = dplyr::recode(EU_2009, "1" = 0, 
                                 "2" = 1, 
                                 "3" = 2, 
                                 "4" = 1, 
                                 "5" = 1, 
                                 "6" = 1), 
         EU_2009 = na_if(EU_2009, 2), 
         Ref_1994 = dplyr::recode(Ref_1994, "1" = 0, 
                                 "2" = 1, 
                                 "3" = 2, 
                                 "4" = 1, 
                                 "5" = 1, 
                                 "6" = 1), 
         Ref_1994 = na_if(Ref_1994, 2),
         Nat_1970 = case_when(Nat_1970 == "1" ~ 0, 
                              Nat_1970 == "2" ~ 1, 
                              Nat_1970 == "3" ~ NA, 
                              Nat_1970 == "4" ~ 1, 
                              Nat_1970 == "5" ~ 1, 
                              Nat_1970 == "6" ~ 1, 
                              Legal_1970 == "1"|Foreign_1970 == "1" ~ NA), 
         Reg_1970 = case_when(Reg_1970 == "1" ~ 0, 
                              Reg_1970 == "2" ~ 1, 
                              Reg_1970 == "3" ~ NA, 
                              Reg_1970 == "4" ~ 1, 
                              Reg_1970 == "5" ~ 1, 
                              Reg_1970 == "6" ~ 1, 
                              Legal_1970 == "1"|Foreign_1970 == "1" ~ NA),
         Mun_1970 = case_when(Mun_1970 == "1" ~ 0, 
                              Mun_1970 == "2" ~ 1, 
                              Mun_1970 == "3" ~ NA, 
                              Mun_1970 == "4" ~ 1, 
                              Mun_1970 == "5" ~ 1, 
                              Mun_1970 == "6" ~ 1, 
                              Legal_1970 == "1"|Foreign_1970 == "1" ~ NA), 
         Nat_1994 = case_when(Nat_1994 == "1" ~ 0, 
                              Nat_1994 == "2" ~ 1, 
                              Nat_1994 == "3" ~ NA, 
                              Nat_1994 == "4" ~ 1, 
                              Nat_1994 == "5" ~ 1, 
                              Nat_1994 == "6" ~ 1), 
         Reg_1994 = case_when(Reg_1994 == "1" ~ 0, 
                              Reg_1994 == "2" ~ 1, 
                              Reg_1994 == "3" ~ NA, 
                              Reg_1994 == "4" ~ 1, 
                              Reg_1994 == "5" ~ 1, 
                              Reg_1994 == "6" ~ 1), 
         Mun_1994 = case_when(Mun_1994 == "1" ~ 0, 
                              Mun_1994 == "2" ~ 1, 
                              Mun_1994 == "3" ~ NA, 
                              Mun_1994 == "4" ~ 1, 
                              Mun_1994 == "5" ~ 1, 
                              Mun_1994 == "6" ~ 1), 
         Nat_2010 = case_when(Nat_2010 == "1" ~ 0, 
                              Nat_2010 == "2" ~ 1, 
                              Nat_2010 == "3" ~ NA, 
                              Nat_2010 == "4" ~ 1, 
                              Nat_2010 == "5" ~ 1, 
                              Nat_2010 == "6" ~ 1, 
                              Eligible_2010 == "3" ~ NA), 
         Reg_2010 = case_when(Reg_2010 == "1" ~ 0, 
                              Reg_2010 == "2" ~ 1, 
                              Reg_2010 == "3" ~ NA, 
                              Reg_2010 == "4" ~ 1, 
                              Reg_2010 == "5" ~ 1, 
                              Reg_2010 == "6" ~ 1, 
                              Eligible_2010 == "2" ~ NA), 
         Mun_2010 = case_when(Mun_2010 == "1" ~ 0, 
                              Mun_2010 == "2" ~ 1, 
                              Mun_2010 == "3" ~ NA, 
                              Mun_2010 == "4" ~ 1, 
                              Mun_2010 == "5" ~ 1, 
                              Mun_2010 == "6" ~ 1, 
                              Eligible_2010 == "2" ~ NA), 
         Nat_2018 = case_when(Nat_2018 == "1" ~ 1, 
                              Nat_2018 == "0" ~ 0, 
                              
                              Eligible_2018 == "3" ~ NA), 
         Reg_2018 = case_when(Reg_2018 == "1" ~ 1, 
                              Reg_2018 == "0" ~ 0, 
                              
                              Eligible_2018 == "2" ~ NA), 
         Mun_2018 = case_when(Mun_2018 == "1" ~ 1, 
                              Mun_2018 == "0" ~ 0, 
                              
                              Eligible_2018 == "2" ~ NA)) %>% 
  mutate(National_total = rowSums(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970")], na.rm = TRUE), 
         National_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970") ])), 
         Regional_total = rowSums(combined[,c("Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970") ], na.rm = TRUE), 
         Regional_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970") ])), 
         Municipal_total = rowSums(combined[,c("Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970") ], na.rm = TRUE), 
         Municipal_eligible = rowSums(!is.na(combined[,c("Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970") ])), 
         EU_total = rowSums(combined[,c("EU_2019", "EU_2009")], na.rm = TRUE), 
         EU_eligible = rowSums(!is.na(combined[,c("EU_2019", "EU_2009") ])), 
         Total_total = rowSums(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                   "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                   "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970", 
                                                   "EU_2019", "EU_2009")], na.rm = TRUE), 
         Total_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                                "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                                "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970", 
                                                                "EU_2019", "EU_2009")])), 
         Domestic_total = rowSums(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                   "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                   "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970")], na.rm = TRUE), 
         Domestic_eligible = rowSums(!is.na(combined[,c("Nat_2018", "Nat_2010", "Nat_1994", "Nat_1970", 
                                                                "Reg_2018", "Reg_2010", "Reg_1994", "Reg_1970", 
                                                                "Mun_2018", "Mun_2010", "Mun_1994", "Mun_1970")])),
                                       National_prop = National_total/National_eligible, 
                                       Regional_prop = Regional_total/Regional_eligible, 
                                       Municipal_prop = Municipal_total/Municipal_eligible, 
                                       EU_prop = EU_total/EU_eligible, 
         Total_prop = Total_total/Total_eligible, 
         Domestic_prop = Domestic_total/Domestic_eligible, 
         Domestic_scaled = scale(Domestic_prop), 
         Total_scaled = scale(Total_prop), 
         Regional_scaled = scale(Regional_prop), 
         Municipal_scaled = scale(Municipal_prop), 
         National_scaled = scale(National_prop), 
         EU_scaled = scale(EU_prop))
                          
#Filtering for complete sets

pairs = combined %>% 
  group_by(LopNrParID) %>% 
  filter(is.na(trust) == F & is.na(behavioral) ==F) %>% 
  filter( n() > 1 ) %>%
  mutate(TWINID = rep(c(1, 2),length.out = n())) %>% 
  ungroup()

#Stacking data by twin pairs

pairs = pairs %>% 
  select(LopNr, LopNrParID, trust, oldtrust, trust_scaled, behavioral, contactpol, contactpub, protest, boycott, finance, petition, TWINID, BESTZYG, SEX, Byear, National_prop, Regional_prop, Municipal_prop, Domestic_prop, Total_prop, interpersonal_1, interpersonal_2, interpersonal_scaled) 

pairs$current_date <- as.numeric(format(Sys.Date(), "%Y%m"))
pairs$current_year <- floor(pairs$current_date/100)
pairs$current_month <- pairs$current_date %% 100

# Calculate age in years and adjust for month of birth
pairs$Byear = as.numeric(pairs$Byear)
pairs$age <- pairs$current_year - floor(pairs$Byear/100) - 1
pairs$age = ifelse(pairs$current_month >= (pairs$Byear %% 100), pairs$age <- pairs$age + 1, pairs$age <- pairs$age)

twin1 = pairs %>% 
  filter(TWINID == 1)

twin2 = pairs %>% 
  filter(TWINID == 2)

final = merge(twin1, twin2, by.x = "LopNrParID", by.y = "LopNrParID", suffixes = c("_T1", "_T2"))

#Creating MZ and DZ datasets
SALTY_MZ <- subset(final, BESTZYG_T1 == 1)
SALTY_DZ <- subset(final, BESTZYG_T1 != 1)

```

```{r Descriptive statistics}

combined %>% 
  filter(BESTZYG == 1) %>% 
 select(trust, oldtrust, SEX, interpersonal_1, interpersonal_2) %>% 
 stargazer(iqr = T, median = T, out = "descriptive_stats_herit_MZ.htm")

combined %>% 
  filter(BESTZYG != 1) %>% 
  select(trust, oldtrust, SEX, interpersonal_1, interpersonal_2) %>% 
  stargazer(iqr = T, median = T, out = "descriptive_stats_herit_DZ.htm")

```

```{r Phenotype Correlations}

pairs %>% 
  select("trust", "behavioral", "contactpol", "contactpub", "protest", "finance", "boycott", "petition", 
         "National_prop", "Regional_prop", "Municipal_prop", "Domestic_prop", "Total_prop") %>% 
  psych::corr.test(ci = T) %>% 
  print(short = F)

total_corr_behav <- cor(pairs$trust, pairs$behavioral, use = "na.or.complete")
total_corr_cpol <- cor(pairs$trust, pairs$contactpol, use = "na.or.complete")
total_corr_cpub <- cor(pairs$trust, pairs$contactpub, use = "na.or.complete")
total_corr_protest <- cor(pairs$trust, pairs$protest, use = "na.or.complete")
total_corr_finance <- cor(pairs$trust, pairs$finance, use = "na.or.complete")
total_corr_boycott <- cor(pairs$trust, pairs$boycott, use = "na.or.complete")
total_corr_petition <- cor(pairs$trust, pairs$petition, use = "na.or.complete")

total_corr_national <- cor(pairs$trust, pairs$National_prop, use = "na.or.complete")
total_corr_regional <- cor(pairs$trust, pairs$Regional_prop, use = "na.or.complete")
total_corr_municipal <- cor(pairs$trust, pairs$Municipal_prop, use = "na.or.complete")
total_corr_domestic <- cor(pairs$trust, pairs$Domestic_prop, use = "na.or.complete")
total_corr_overall <- cor(pairs$trust, pairs$Total_prop, use = "na.or.complete")

```

```{r Phenotype correlations table}

# create empty data frame to store results
corr_df <- data.frame(
  "Trust" = character(),
  "Variable" = character(),
  "Correlation" = numeric(),
  "95% CI" = character(),
  stringsAsFactors = FALSE
)

# iterate over variables

vars <- c("behavioral", "contactpol", "contactpub", "protest", "finance", "boycott", "petition", "National_prop", "Regional_prop", "Municipal_prop", "Domestic_prop", "Total_prop")

for (var in vars) {
  
  # run correlation test
  corr_test <- cor.test(pairs$trust, pairs[[var]], method = "pearson", 
                         use = "pairwise.complete.obs")
  
  # extract values
  cor_coef <- round(corr_test$estimate, digits = 2)
  ci_lower <- corr_test$conf.int[1]
  ci_upper <- corr_test$conf.int[2]
  
  # add results to data frame
  corr_df <- rbind(corr_df, data.frame(
    "Variable" = var,
    "Correlation" = cor_coef,
    "95% CI" = paste0("(", sprintf("%.2f", ci_lower), ", ", sprintf("%.2f", ci_upper), ")")
  ))
}

# print table
rownames(corr_df) = c(1:12)
print(corr_df)

corr_df = corr_df %>% 
  pivot_wider(names_from = Variable, values_from =  c(Correlation, X95..CI))

```

```{r Political Trust ACE, echo = FALSE}

#Political trust
ACE <- umxACE("ACE", selDVs = "trust_scaled", selCovs = c("SEX", "age"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, optimizer = "SLSQP", tryHard = "yes")
set.seed(1234)
ACE <- umxConfint(ACE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(ACE, verbose = TRUE)

AE = umxModify(ACE, update = c("c_r1c1"), comparison = TRUE)
AE <- umxConfint(AE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(AE, verbose = TRUE)

CE = umxModify(ACE, update = c("a_r1c1"), comparison = TRUE)
CE <- umxConfint(CE, run = TRUE, level = .95, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(CE, verbose = TRUE)

umxCompare(ACE, comparison = c(AE, CE))

```

```{r Interpersonal Trust ACE, echo = FALSE}

#Political trust
ACE <- umxACE("ACE", selDVs = "interpersonal_scaled", selCovs = c("SEX", "age"), sep = "_T", mzData = SALTY_MZ, dzData = SALTY_DZ, optimizer = "SLSQP", tryHard = "no")
set.seed(1234)
ACE <- umxConfint(ACE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(ACE, verbose = TRUE)

AE = umxModify(ACE, update = c("c_r1c1"), comparison = TRUE)
AE <- umxConfint(AE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(AE, verbose = TRUE)

CE = umxModify(ACE, update = c("a_r1c1"), comparison = TRUE)
CE <- umxConfint(CE, run = TRUE, level = .95, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(CE, verbose = TRUE)

umxCompare(ACE, comparison = c(AE, CE))

```

